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This paper addresses the question of how Brownian-like motion can arise from the solution of 
a deterministic differential delay equation. To study this we analytically study the bifurcation 
properties of an apparently simple differential delay equation and then numerically investigate the 
probabilistic properties of chaotic solutions of the same equation. Our results show that solutions of 
the deterministic equation with randomly selected initial conditions display a Gaussian-like density 
for long time, but the densities are supported on an interval of finite measure. Using these chaotic 
solutions as velocities, we are able to produce Brownian-like motions, which show statistical prop- 
erties akin to those of a classical Brownian motion over both short and long time scales. Several 
conjectures are formulated for the probabilistic properties of the solution of the differential delay 
equation. Numerical studies suggest that these conjectures could be "universal" for similar types of 
"chaotic" dynamics, but we have been unable to prove this. 



I. INTRODUCTION 

In 1828, Robert Brown reported his observations of the 
apparently erratic and unpredictable movement of small 
particles suspended in water, a phenomena now known as 
"Brownian motion" . Almost three-quarters of a century 
later, a theoretical (and essentially molecular) explana- 
tion of this macroscopic motion was given by Einstein, in 
which Brownian motion is attributed to the summated 
effect of a vary large number of tiny impulsive forces de- 
livered to the macroscopic particle being observed [l| (A 
nice English translation of this, and other works of Ein- 
stein on Brownian motion can be found in Fiirth [2[). 
Brownian motion has played a central role in the model- 
ing of many random behaviors in nature and in stochas- 
tic analysis, and formed the basis for the development of 
an enormous branch of mathematics centered around the 
theory of Wiener processes. 

Since Brownian motion is typically explained as the 
summated effect of many tiny random impulsive forces, 
it is of interest to know if and when Brownian motion can 
be produced from a deterministic process (also termed as 
deterministic Brownian motion) without introducing the 
assumptions typically associated with the theory of ran- 
dom processes. Studies starting from this premise have 
been published in the past several decades, and there 
are numerous investigations that have documented the 
existence of Brownian-like motion from deterministic dy- 
namics, both in discrete time maps and flows . These 
models have included the motion of a particle subjected 
to a deterministic but chaotic force (also known as micro- 
scopic chaos) 0, m, or a many-degree-of-frcedom Hamil- 
tonian Experimental evidence for deterministic 

microscopic chaos was reported in Q by the observation 
of Brownian motion of a colloidal particle suspended in 
water (cf. [l3| for a more tempered interpretation, and 



[111 Chapter 18] for other possible interpretations of ex- 
periments like these). 

Several investigators have shown that a Brownian-like 
motion can arise when a particle is subjected to impul- 
sive kicks, whose dynamics are modeled by the following 
equations d [H [II 



dx 

lit 
dv 
I— 
dt 



~iv + /(<)• 



(1) 



In equation / is taken to be a fluctuating "force" 
consisting of a sequence of delta-function like impulses 
given by, for example. 



/(t) = mK^C(t)5(i~nr), 



(2) 



n=0 



and ^ is a "highly chaotic" deterministic variable gener- 
ated by ^(^ -I- r) = r(^(t)), where T is an exact map or 
semi-dynamical system, e.g. the tent map on [—1, 1] (for 
more discussions and terminologies see 0, and refer- 
ences therein). In the equations the impulsive 
forces are described by £^{t)5{t — nr), which are assumed 
to be instantaneously effective and independent of the 
velocity v{t). Dynamical systems of the form ([T]) have 
received extensive attention, and are known to be able 
to generate a Gaussian diffusion process H, [l^-EBl • 

In this study, we sought an alternative continuous time 
description of the "random force" f{t), which was as- 
sumed to depend on the state (velocity) of a particle, 
but with a lag time r, i.e., 

fit)^F{vit-T)), (3) 

and where F has the appropriate properties to generate 
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chaotic solutions. Thus, we consider the foUowing differ- 
ential delay equation 



dx 



n— = -jv + F{v{t - t)), 
v{t) = (/)(t), -r < i < 0, 



(4) 



where (f)(t) denotes the initial function which must al- 
ways be specified for a differential delay equation. The 
second equation in (|4|) is known to have chaotic solutions 
for some forms of the nonlinear function F, for exam- 
ple see [l6l - [20j . In these cases, the force F{v{t — t)) 
is certainly deterministic but is also unpredictable (in 
practice, but not in principle) given knowledge of the 
initial function. In this paper, we will examine how a 
Brownian motion can be produced by the differential delay 
equation In particular, we investigate the statistical 
properties of the velocity v{t), and ([Aa;(t)]^), the mean 
square displacement (MSD), of the solutions defined by 
dH). We note that unlike the equation ([T|), which is lin- 
ear and non-autonomous, the equation (0]) is a nonlinear 
autonomous system. Numerical simulations have shown 
that the second equation in @ can generate processes 
with a Gaussian-like distribution [l^, [T^. Nevertheless, 
to the best of our knowledge, there is no analytic proof 
for the existence of Brownian motion based on the differ- 
ential delay equation 

We first make some observations about the second 
equation in (|4]) which determines the dynamics of the 
velocity. A simple form of the "random" force is binary 
and fluctuates between ±/o, for instance, given by 



F{v) = 2/o 



H{sin{2TTfiv)) - \ 



where H is the Heavyside step function, i.e., 

„, . f for w < 
= I 1 for « > 0. 

Then we have following equation 
dv 



dt 
v{t) 



-jv + 2 {H{sm{2nl3v{t - 1))) 
(f>it), -T <t<0. 



(5) 



(6) 



(7) 



Here and later we always assume the mass m = 1 and 
/o = 1 which can be achieved through the appropriate 
scaling. The delay differential equation ([7]) with a binary 
"random force" can be solved iteratively by the method of 
steps [2l| . Despite its simplicity, it can display behaviors 
similar to a random process. An example solution of ([7]) 
is shown in Figure [TJ which looks like noisy. 

The "random force" in ([7]) is discontinuous and gives 
a continuous zigzag velocity curve (c.f. inset in Figure 
[T|). In this paper, we will instead study an analogous 
different differential delay equation 



dv 

Itt ~ 
v{t) 



— 7W - 



s\n{2TTl3v{t 
-\<t<Q. 



1)), 



(8) 
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FIG. 1. (Color online) A sample solution of Q with /3 = 
10, 7 = 1, and an initial function 4>{t) = —0.1, t £ [—1, 0]. The 
rectangular inset shows the solution segment for 98 < t < 100. 



In (|8]), the parameter /3 measures the frequency of the 
nonlinear function, and will turn out to be an essential 
parameter in the present study. Note that one can re- 
scale and translate the variables such that equation ([5]) 
can be rewritten as 



dv ■ r r \ \ 

— = -V + fism(v[t - r) - xq), 
dt 

v{t) = <f){t), -r < t < 0. 



(9) 



Equation ([9|) (also known as Ikcda equation) was pro- 
posed by Ikeda et. al. to model a passive optical bistable 
resonator system, and shows chaotic behaviors at partic- 
ular parameters such as /.( = 20, = 7r/4 and t = 5 

Hi!. 

In this paper, we will study the dynamical properties 
of the solutions of ([8]), both analytically and numerically. 
We focus in particular on the probabilistic properties of 
the chaotic solutions. Then we investigate chaotic solu- 
tions of 



dx 

— ~ v 
dt 

= -jv + sm{2TT/3v(t~l)), 
v{t) = 4>{t), -1 < i < 0, 



(10) 



and characterizing the statistical properties as com- 
pletely as we can. The main result is to show that 
the equation (|10p can reproduce experimentally observed 
data of Brownian motion over a wide range of time scales, 
in spite of the fact that the evolution equation is deter- 
ministic. Therefore, deterministic Brownian motion can 
be generated from the equation (fTOj) . 

The outline of the rest of this paper is as follows. We 
first perform a bifurcation analysis for equation ([5]) in 
Sectionini In Section lllll we study the probabilistic prop- 
erties of the chaotic solutions numerically. In Section HVl 
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we numerically examine the dynamics of the chaotic so- 
lutions of (jlOp and compare our results with recent ex- 
perimental measurements of the motion of a Brownian 
particle [2^. Section fVl presents five conjectures based 
on our studies that we have been unable to prove, but 
which we believe to be true. These conjectures indicate 
a possible direction for the analytical proof of the exis- 
tence of deterministic Brownian motion from differential 
delay equation Finally, we conclude the paper with 
discussion and conclusions in Section IVll 

II. BIFURCATION ANALYSIS 

In this section, we commence our study by performing 
a bifurcation analysis for equation ([5]). Wc always assume 

7 = 1- ^ ^ 

The bifurcation structure of the Ikeda equation has 
been studied several times from different perspectives 
[23 - |29j . Here, we present a complete picture (see The- 
orem [T] below) for the bifurcation structure of the equa- 
tion dS]), which has not, to the best of our knowledge, 
appeared previously. 

A. Steady state solutions 



to. In this case, {/3,v*) is a Hopf bifurcation point of ([5]). 
(Throughout this paper we will use the notation (/3,w*) 
where /3 is the bifurcation parameter and v* is the bifur- 
cation point.) Figure [2] graphically displays the steady 
states for < /3 < 6. 
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FIG. 2. (Color online) Steady states v* of ^ for different 
values of /3 (here 7=1). The inset shows the detail in the 
small rectangular area. Black denotes locally stable steady 
states (ss), while green (light gray) denotes unstable steady 
states (us). 



The steady states of equation ^ are given by the so- 
lutions of 



V = sin(27r/3z;). 



When /3 < l/(27r), equation ([TT|) has only one real solu- 
tion, namely v = 0. When /3 > l/(27r), ([TT]) has (4[/3]-hl) 
real solutions where [/3] denotes the integer part of /?. 
These solutions are separated by critical points which 
are given by the roots of 

1 = 27r/3cos(27r^u), 



I.e. 



k 1 f 1 

t: ± 7^ — t: arccos - — - 
l3 27r/3 \2nl3 



<k < 



Let V* be a steady state of Linearization of ([5]) 
around v ~ v* gives 

dv 

— = -■D + 27r^cos(27r^w*)S(t- 1). (12) 
at 

Thus, the steady state solution v{t) = v* is locally stable 
if and only if 



secw < 27r^cos(27r^u*) < 1, 
where lo + tancj = 0, cj G (0, tt). 



(13) 



In particular, if 27r/3 cos(27r/3ti*) = secw, the linearized 
equation (jl2[) has a pair of complex conjugate eigenval- 
ues, and therefore has a periodic solution with frequency 



We will first state the following Lemma before giving 
the main results of our bifurcation analysis. 



(11) Lemma 1 For any a £ (— oo, 0) U (^, -foo), let 



p(«^2,Jl-(^) -«.■»« (^), (14) 



Then for any fc G N* , the equation 
27rfc = piP) 



(15) 



has a unique solution /? > |a|/27r. In particular, when 
a — 1 and k — 0, we have /? — l/27r. 



Proof. Since P^-^) < 0, and when P > ^ 



47r/32 - a a{2na - 1 



> 



>0, 



the first part of the Lemma follows. It is easy to verify 
that when a = 1, then p(^) = 0. ■ 

In the following, we define sequences {ofe} and {6fe} 
such that 



2-Kk = 2-Kak\l 1 ~ ( ^ I - arccos ( ^ 1 , fc e N*, 



V 2'Kai 



2'Kai. 



(16) 



and 



27rfc = 2TrbkA 1- 



2TTbk 



secw \ / secw , 

— arccos | „ . | , KG 



27r6t 



(17) 
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Note that oq = l/27r. Furthermore, since sec^w > 1 and 
arccos(^^) > arccos( j;^^), we have 



^ n L. ( sccwX / scca;\ 

and therefore 6^ > a^. When fc — > oo, the solution of 
([T5|) is approximately 

1 afa — 2) 



Thus, for large fc wc have 

(sccw — 1)^ 



- flfc 



47r2fc 



(18) 



Note that the sequences {au\ and {6^} can be ordered as 

< flo < &o < fli < &i • • • < Ofe < 6fc < • • • . (19) 

In the following Theorem, we prove that {a/c} are saddle 
node bifurcation points and {fofc} are Hopf bifurcation 
points of Eq. 

Theorem 1 Consider Eq. (|5]) and its steady state solu- 
tions. Let {ak} and {bk} be defined as above. Then: 

(1) When < /? < flo, (O has only one steady state 
V* = 0, and it is locally .stable. 

(2) When j3 = oq, {P,v) = (ao7 0) is a pitchfork bifur- 
cation point of (O . 

(3) When (3 > oq, the steady state v* = is unstable. 

(4) For any k eW, let 



1 



yk = — 



Ok 2nak 



27rai 



(20) 



Then (13, v) — {ok^^yk) are saddle node bifurcation 
points of 



(5) For any k eW, let 



k 1 ( secoj 

bk 27r6fc V 27r6fe 



(21) 



Then (/3,w) = {bk,izk) are Hopf bifurcation points 
of ®. 

(6) For every k € Nq, and ak as defined above, there 
is a function fk{P) which is continuous on [0^,00), 
such that fk{ak) = yk, o.nd when /3 > Ok, v ~ 
ztfkiP) satisfies ([TT]) . and 



l-2^/3cos(27r/3/fc(/3)) <0. 



(22) 



The steady state solutions v{t) = ±fk{f3) are un- 
stable. 



(7) For every fc € Nq, and ak, bk as defined above, there 
is a function gk{P) , which is continuous on [ak,oo), 
such that gk{ak) = yk, gk{bk) = Zfc, and when /3 > 
ak, V = ±gk{l3) satisfy ([TT|) . and 



1 - 27r/3cos(27r/35fc(/3)) > 0. 



(23) 



Further, 



(a) When ak < (3 < bk, the steady state solutions 
v{t) = ±gk{P) are locally stable. 

(b) When (3 > bk, the steady state solutions v{t) = 
±gk{f3) are unstable. 

(8) When (3 increases past the Hopf bifurcation point 
bk, the two steady state solutions v{t) = ±gk{l3) 
lose stability and generate a periodic solution, with 
angular frequency lo. Therefore there exists a se- 
quence {ck}, such that the periodic solutions gen- 
erated from the Hopf bifurcations are stable when 
bk < 13 < Ck. 

Proof of Theorem [T] is given in the Appendix. 

From Theorem [U when /3 increases from 0, in any in- 
terval /? G {ak,bk), equation ([5]) has two stable steady 
state solutions. The length of the intervals (a^, bk) tends 
to zero as 1/fc as fc — > 00. In other situations, however, 
all steady states are unstable, and therefore complicated 
dynamical behaviors may be expected. An exploration 
of the nature of these constitutes the remainder of this 
paper. 



B. Periodic solutions and chaotic attractors 

In this section, we numerically investigate the long 
term behavior of the solutions of equation ([5]). For each 
value of the parameter /? > 0, we solve equation ^ to 
obtain 100 independent sample solutions, each with ran- 
domly selected constant initial function 

vit)^vo {-l<t<0), 

where vq € (—1,1) and is uniformly distributed. Each 
solution is obtained using Euler's method (with a time 
step At = 0.001) up to t = 500, such that the solution 
reaches a stable state (either oscillatory or steady state) , 
and the resulting data from 300 < t < 500 are used for 
further analysis as detailed below. 

To distinguish oscillatory solutions from constant so- 
lutions in the simulation, wc investigated the upper and 
lower bounds of v{t) in 300 < t < 500, denoted by Umax 
and Umin, respectively. Therefore, a solution is considered 
to have approached a stable steady state if Wmax — '^min, 
and approached a stable oscillatory solution between Vmin 

and Umax if Wmin < Wmax- 

Figure [3] shows the simulation results. For each value 
of /3 € (0,3), there are 100 pairs of dots from 100 ini- 
tial functions, corresponding respectively to Umin (red 



(light gray) dots) and Vmax (blue dots) of a solution. 
Stable steady state solutions are shown by the super- 
position of blue and red (light gray) dots. We are inter- 
ested in the oscillatory solutions, which are indicated by 
well separated blue and red (light gray) dots in Figure [31 
There are two types of oscillatory solutions. Regular pe- 
riodic solutions appear when /? is close to (but greater 
than) the Hopf bifurcation points bk- The amplitude 
of these regular solutions only depend on the parame- 
ter /3. Irregular oscillatory solutions occur for almost all 
/3 (/3 > 0.85 except small gaps at /3 G (1.019, 1.033) and 
/? G (1.270, 1.360), respectively. The amplitude of these 
irregular oscillatory solutions depend on the parameter 
/3 as well as the initial functions. 

The bifurcation diagram in Figure [3] shows clear evi- 
dence for multi-stability of the solutions of equation ([5]) : 

1. When ao < (3 < bo, there are two stable steady 
states. 

2. When bo < 13 < Co, there are two stable periodic 
solutions. 

3. When ak < (3 < bk{k G N*), there are two stable 
state states, and stable irregular oscillations. 

4. When bk < P < Ck{k G N*) there are two stable 
periodic solutions, and stable irregular oscillations. 

Some sample solution examples for different values of (3 
are shown in Figure HJ and these also illustrate the ex- 
istence of multi-stability of solutions dependent on the 
initial function. 



III. PROBABILISTIC PROPERTIES 

In the numerical bifurcation analysis of Section [III we 
have shown that when 7 = 1 and (3 > 0.85, the differen- 
tial delay equation (|8]) has irregular oscillatory solutions 
that display chaotic behavior. In this section, we numer- 
ically study the probabilistic properties of these irregular 
solutions. We will show in Section fill Fl that when f3 is 
sufficiently large, on a sufficiently long time scale these 
chaotic solutions behave like a noise source with a "trun- 
cated" Gaussian density. 

A. Numerical scheme 

Throughout this section, the probabilistic properties of 
solutions of equation ([5]) are studied numerically. In the 
numerical simulations, for a given set of parameters, we 
solve the equation ([5]) with a randomly selected constant 
initial function 

v{t)=voe{-l,l), i-l<t<0), (24) 

where vo is drawn from a uniformly distributed density. 
The solution v{t) is solved using Eulcr's method (with 
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FIG. 3. (Color online) Bifurcation diagram of ([H]) for different 
/3 (here 7 = 1). Here PF denotes the occurrence of a pitch- 
fork bifurcation, HP stands for a Hopf bifurcation, and SN a 
saddle node bifurcation. In the diagram, black dots denote 
stable steady states (ss), green denotes unstable steady states 
(us) (also refer to Figure [2]), blue is used to indicate the up- 
per bound (ub) of oscillatory solutions, and red (light gray) 
indicates the lower bound (lb) of oscillatory solutions. 



a time step At ~ 0.001) up to t = 10^, and is sampled 
every lO'^ steps to generate a time series where 
Vn = v{n X lO'^At). The resulting time series of values 
{«„} is used to characterize the statistical properties of 
the solution. In particular, we focus on the mean value 
/J,, the upper bound K, the standard deviation a, and 
the excess kurtosis 72 of the time series, which are re- 
spectively defined by 

1 ^ 1 ^ 

n—1 n—1 

1 " 

72 = ^ - 3, where = 77 y](t'« - m)''- 

n—1 

(25) 

The excess kurtosis 72 measures the sharpness of the den- 
sity of the sequence, and a value of 72 = is characteristic 
of a normal Gaussian distribution. 

In the following discussion, we will show that when (3 
is outside the region of bistability, i.e., 

/3€/= (^,+oo]\[J[afc,Cfc], (26) 
^ ^ fc=0 

the statistical properties are independent of the initial 
function vo and sampling frequencies. Therefore, the 
quantities defined by (j25p only depend on the parame- 
ters 13 and 7, and this dependence is discussed below. 
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p=1.2 p=1.3 p=1.45 P=1.75 




J ^ [ J ^ [AAAAAAAA/^ ^ [ . . . 

10 20 30 10 20 30 10 20 30 10 20 30 
t t t t 



FIG. 4. (Color online) Sample solution segment examples for different values of /3 (from small to large) as indicated in each 
panel. The time point in figures correspond to t — 400 in the simulations. Each panel contains between one and four solution 
curves (marked by different colors and also locations), each of which corresponds to a stable oscillatory solution. For a given 
values of /?, the initial functions of these solutions are different from each other, and taken as constant in the interval — 1 < t < 0. 
The constant initial functions are: /3 = 0.6: vq = 0.55, -0.55; /3 = 0.77: -yo = 0.41129, -0.41129, 0.630575, -0.630575; ^ = 1.0: 
vo = 0.55; P = 1.025: vo = 0.55; /3 = 1.2: vo = 0.55; j3 = 1.3: vo = 0.55; /3 = 1.45: vo = 0.975, 0.1, -0.975; /3 = 1.75: vo = 0.1. 
Here 7 = 1. 



B. Stationary density of solutions 

From the bifurcation analysis in Section |IT1 equation 
^ displays bistability when /? G [a/j,Cfc] for some A: e N. 
This suggests that when 

P eh=[j[ak,CK], (27) 

km 

a solution of ([5]) with a randomly selected initial func- 
tion will converge to one of the stable branches (either a 
steady state, or a periodic solution, or an irregular oscil- 
latory solution) as shown by Figure |31 Thus, the station- 
ary density of all solutions of ^ is expected to be multi- 
modal, and therefore the limiting statistical properties of 
a solution (as mirrored in the density constructed along 
the solution trajectory) depend on the initial condition. 
The upper two panels of Figure [5^ show the multi-modal 
distributions of the stationary densities when /3 = 1.25 
and /3 = 1.35, respectively. The results are obtained from 
10^ independent solutions at t = 100, each with a ran- 



domly selected constant initial function as given by (|24p . 

In the following discussion, we focus on the alternative 
situation in which 

/3g/= (^^,c»^\/o. (28) 

In this case, the numerical bifurcation analysis shown in 
Figure [3] indicates that any solution of ([8]) converges to 
an irregular oscillatory solution irrespective of the initial 
function. (Note that solutions with different initial func- 
tions will not, in general, converge to the same solution 
as can be seen by multiple values of the upper and lower 
bounds of the oscillatory solutions in Figure [3]) Despite 
the fact that the solutions in such situations are not the 
same, we will see that these solution trajectories share 
the same statistical properties. Figure [5^ (black curve in 
the lowest panel) shows the stationary density obtained 
from 10^ independent solutions (here /3 = 2.0 £ /), each 
with a randomly selected constant initial function (j24p . 
The result is a uni-modal density of the distribution of 
solution values along the trajectory. Alternatively, the 
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same density can also be obtained through a time series 
{vn} of a solution with randomly selected constant initial 
function (Figure [5^, red (light gray) curve in the lowest 
panel, refer to Section [III Al for details). Furthermore, 
we have found that the same density of the distribution 
is obtained when we choose different forms for the initial 
function, such as a sinusoidal initial function, or a poly- 
nomial initial function (data not shown). These results 
strongly suggest that these irregular solutions are ergodic 
in some sense, i.e., the statistical properties of one solu- 
tion are the same as those of an ensemble of independent 
solutions [30j . 

In what follows, we focus on quantifying the statistical 
properties of these irregular solutions. An example of one 
of these irregular solutions is shown in Figure [SJd, with 
the corresponding power spectrum w(/) shown in Figure 
[S];. Note that the power spectrum is essentially flat with 
no predominant characteristic frequency, indicating that 
the solution is, indeed, chaotic. 



C. Dependence of the statistical properties on /3 

Now, we take 7 = 1, and /3 £ [1,50] to numerically 
study the statistical properties of the irregular solutions. 
For each /3 the sampled time series of a solution v{t) 
is used to obtain the mean value fi, the upper bound K, 
the standard deviation cr, and the excess kurtosis 72 (refer 
Eq. ([^5]) for these definitions). Figure [5] shows these four 
statistical indicators as functions of the parameter (3. 

Figure [6^ shows the mean value as a function of /3, 
indicating that /i(/3) ~ 0. Figure [BJd shows the bound K 
as a function of /3. The numerical results show that K 
decreases with /3, and can be accurately approximated 
by 



0.02 



K 



1 



0.68^^+0.60' 



(29) 



Note that as /3 — > 00. solutions of ^ are bounded and K 
varies as 

/3-1/2. Fi 

gureint shows the standard deviation 
cr as a function of (3. The standard deviation decreases 
with /3, and can be fitted with 



0.32/3 



'1/2 



(30) 



As in the case of the upper bound, when /3 00, we also 
find that a varies as f3~-^/^. 

Figure |6ji shows the excess kurtosis 72 as a function 
of /3, and the numerical results reveal that the excess 
kurtosis increases with /3 towards 0, approximately as 
— 1/(3. The negative value indicates that the distribu- 
tion is platykurtic (the tail of the distribution is thinner 
relative to a Gaussian). This is because the solution is 
bounded, and therefore the tail is truncated. Note that 
a larger /? means a smaller absolute value of the excess 
kurtosis, and thus that the distribution is more like a 
Gaussian distribution. 
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FIG. 5. (Color online) (a) Stationary density functions of 
solutions of (O for different values of /3 (as shown in the pan- 
els). The density functions are obtained from 10* independent 
solutions, each with a randomly selected constant initial func- 
tion as in Eq. ()24|) . The red (light gray) curve in the lowest 
panel is the density function obtained from a random solution 
trajectory with a constant initial function, (b) Example of a 
segment of a single solution with /3 = 2.0, and initial function 
vq — 0.04. (c). Power spectrum of the solution in (b). In the 
simulations, 7 = 1.0. 



D. Dependence of the statistical properties on 7 

We now fix /3 = 20 and study the dependence of the 
statistical properties on 7. Figure [7] shows the simulation 
results with the same statistical indicators plotted as in 
Figure [6] The red (light gray) curves in Figure [7] are fit 
by 



A-(/3,7) 
72 (/3, 7) 



1 



V7(0.68V? + 0.60^7) ^ 
0.32 
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(31) 

(32) 
(33) 
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FIG. 6. (Color online) (a) Mean value /i as a function of /3. (b) 
The upper bound K as a function of /3. The solid curve shows 
the fit to iC = 1/(0.687^ + 0.60). (c) Standard deviation a as 
a function /3. The solid curve is the graph of a = 0.32/3"^'^^. 
(d) Excess kurtosis 72 as a function of /3. Solid curve shows 
the fit to 72 — — 1//3. Remember that for all of these results, 
7 = 1, and the initial functions are constants as in (|24p . 



FIG. 7. (Color online) (a) Mean value as a function of 7. 
(b) The upper bound K as a function of 7. (c) Standard 
deviation cr as a function 7. (d) Excess kurtosis as a function 
of 7. Solid curves in (b)-(d) show the fits (|31l) - (|33| ). Here 
13 = 20 and the initial functions are constants as in Eq. (|24|) . 
Results for 7 > 30 are not shown due to numerical instability. 



The functions (PT|) - (P5|) give the general dependence of 
the statistical indicators with equation parameters /3 and 
7, and are obtained as follows. First, we rescale equation 
([8]) by introducing u — 7w,/3' = [3/ J- Then u{t) satisfies 



du 
Hi 



= 7(-u + sin(27r^'u(i - 1))). 



(34) 



The statistical indicators K , a and 72 for solutions u{t) of 
P4p are independent of 7 (data not shown) , and depend 
on P' through the same functions as in Figure El There- 
fore, we obtain the functions (|3T|) - (|33)) with the scaling 
V = M/7. 



E. Correlation function 

Next, we investigate the correlation function of a solu- 
tion of the differential delay equation (jS]). The normal- 
ized correlation function of a solution is defined as 



C(r) 



lim 



v{t)v{t + r)dt 



(35) 



Figure [8ji shows the correlation function C(r) for dif- 
ferent values of (3 (with 7 = 1). From Figure |S1 the cor- 



relation function can be approximated as an exponential 
function of the form 



C{r) 



-r/to 



where the constant to gives the correlation time. Figures 
[5]3-c show that the correlation time is largely independent 
of /?, and that it is approximately given by I/7. 



F. Quasi-Gaussian distribution 



From our numerical results in Sections IIII Cl and llll Dl 
it is clear that the excess kurtosis 72 of one of the irregular 
solutions of ([8]) varies with /? and 7 according to 72 — 
— 7//?. Thus, the distribution approaches a Gaussian-like 
distribution when /3 is large (and 7 is fixed), but one with 
a truncated tail so that it is supported on a set of finite 
measure. We call such a truncated Gaussian distribution 
a quasi- Gaussian distribution, and consider these further 
in this section. 

Let /i and a be the mean and standard deviation of a 
quasi-Gaussian noise, and assume that the noise signal is 
supported on an interval [fi — K,fi + K]. Then the density 
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FIG. 8. (Color online) (a) Correlation function C(r). Here, 7=1, and /3 = 5 (blue circles), 10 (red up triangles), 15 (black 
down triangles), 20 (green squares), respectively, (b) Correlation time as a function of /3 (with 7 = 1). (c) Correlation time as 
a function of 7 (with /3 = 20), solid curve is the fit with to = 1/7- 




1 



27ra{'l>{K/a) - ^{-K/a)) 



if \v-n\<K (-36) 
other wise, 



(37) 



and 



1 



(38) 



In particular, when fi = and cr = 1, we have a standard 
quasi-Gaussian distribution, with density function 



p{v;0,l,Ko) 



0, 



\v\ < Ko 
other wise. 



where 



Co 



J —Kn 



(39) 



(40) 



There is only one adjustable parameter, namely the 
bound Kq, in a standard quasi-Gaussian distribution. 

We can now compare the distribution function ob- 
tained from our simulation data with the quasi-Gaussian 
distribution. To this end, we first normalized the signal 
sequence {«„}. From Section UlI Dl let = i;„/ct(/3,7) 
so the sequence {Cn} has mean fi = 0, standard deviation 
CT = 1, and is bounded by 



^(/3,7) 



(41) 



Equation (|4T|) gives the relation between the equation 
parameter /3/7 and the adjustable distribution parameter 



Kq, which is shown in Figure [Hti. We have Kq ~ 5 when 
/3/7 is large. 

Figure shows the result of fitting the density func- 
tion Eq. ([55]) with our simulation data. We can see 
that (123) provides a reasonable fit for the simulation 
data when /3 > 6 (and 7 = 1). Thus, the simulation 
results indicate that when /3/'-f is large, the density of 
the distribution of irregular solution trajectories of the 
differential delay equation ([5]) can be approximated by a 
quasi-Gaussian distribution. 



IV. DETERMINISTIC BROWNIAN MOTION 

Now, we will show that the differential delay equation 
pop , in a suitable parameter region, can generate dy- 
namics with many of the properties of Brownian motion 
in spite of the fact that the evolution equation is de- 
terministic. Therefore, these dynamics are examples of 
deterministic Brownian motion. 

Gonsider solutions of the following deterministic sys- 
tem: 



dx 
dt 

= -7W + sin(27r/3w(t - 1)) 



(42) 



Here x(t) measures the position of a particle with veloc- 
ity v{t). Figure [TUk shows sample solutions x{t) of ([^ . 
which are akin to the dynamics of a Brownian particle. 
Figure [TOb shows the density functions of Ax{t), the dis- 
placement from the particle initial position at different 
times t. These numerical results show that at any time 
t, Ax{t) has Gaussian like distribution. 

Figure [TTI shows {[Ax{t)]'^), the dependence of mean 
square displacement (MSD), as a function of t. In the 
simulations, we set B = 20. and chose different values 
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FIG. 9. (Color online) (a) The normalized bound K as 
a function of /3/7. The solid curve shows the fit with 
\//3/7/(0. 21-1^/3/7 + 0.19). Circles are data from simulations 
in Sections IIII Cl and llllDI (b) Density functions obtain from 
a solution of the differential delay equation with given value 
of P (and with 7 = 1). Solid curve shows the density function 
of quasi-Gaussian distribution according to Eq. H39[) and with 
Ko obtained from j3 = 10. 



of 7 (as shown in the figure panel). For each 7, the 
MSD ([Aa;(t)]^) is obtained from 10'^ independent tra- 
jectories, each with a randomly selected constant initial 
velocity. In Figure [TTl we normalized the results for dif- 
ferent parameters through D = cr'^ /j, the "diffusion con- 
stant". Our simulations show that {[Ax{t)]'^) = 2Dt at 
long time scales, as predicted by Einstein't theory for 
Brownian motion [l|. At short time scales, we have 
([Aa;(t)]^) ~ cDt^, where the pre-factor c depends on 
the initial condition. This result agrees well with re- 
cent measurements of Brownian motion using an optical 
tweezer [2^. In the interpretation of experimental data, 
it is a long debated question if and by what mens we can 
distinguish whether an observed irregular signal is de- 
terministically chaotic or stochastic [3l|, |32| . Results in 
the current study indicate that experimentally observed 



FIG. 10. (Color online) (a) The numerically produced deter- 
ministic Brownian motion x{t). Here /3 — 20,7 = 1- Initial 
velocities are constants as in Eq. (|24|) . (b) Density func- 
tions of Ax{t) (= x{t) — x{0)) at different times t. The sym- 
bols are taken from the numerical solutions, while the solid 
curves show the density function of the corresponding quasi- 
Gaussian distribution. 



data for Brownian motion can be reproduced by solutions 
of a deterministic differential delay equation over a wide 
range of time scales (six orders of magnitude). Thus, dif- 
ferential delay equations provide an alternative way for 
reproducing "random" signals. 



V. CONJECTURES 

The results that we have presented to this point are 
so intriguing that we are led to formulate a series of con- 
jectures. Though we believe these to be true, all efforts 
to prove them have proved fruitless to date. We present 
them in the hope that others will find their proof a chal- 
lenge that they are able to overcome. 

In formulating these conjectures, we focus on the irreg- 
ular solutions for large (3, and therefore we will always 
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FIG. 11. (Color online) Mean square displacement (MSD) 
{[Ax{t)]^) of deterministic Brownian motions. Here D = 
which is analogous to the diffusion constant of a 
Brownian particle in solution. The dashed lines show 
{[Ax{t)]'^)/D = at short times, and {[Ax{t)]^) = 2t for a 
longer time scale respectively. Here /3 = 20 and 7 are shown 
in the figure panel. 



Remark 1 When 13 ^ I, we have to exclude the cases in 
which Vfj{t; (j)) converges to either a stable steady state or 
a stable periodic solution arising through a Hopf bifurca- 
tion. Thus, with this exclusion, the solution vi3{t] </>) will 
converge to either a periodic solution which is symmetric 
about 0, or an irregular solution. Equation (09]) always 
holds for a symmetric solution. Therefore, to prove Con- 
jecture [T] wc only need to show that (HHJ is satisfied by 
the irregular solutions for any /3 > 0. 

Conjecture 2 For any (j) G C([-1,0],M) satisfying 
the limit 



lim P^''^K{p- (j)) 

l3eI,l3-foo 

exists, independent of (j), and is positive. 



(50) 



Conjecture 3 For any (I) & C([— 1,0],M) satisfying pS)) . 
the limit 



lim (j)) 

0^1,13-^00 

exists, independent of (j), and is positive. 



(51) 



assume that (3 is such that equation ^ has no stable 
steady state or stable periodic solution. In particular, 
according to Theorem [U we will always assume 7 = 1 
and /3 el, with / defined by (pS]) . 
Let Vij{t; (j>) be the solution of 



dzj 

= -v + sm(2nl3v(t - 1)), 

at 



-l<t<0. 



Define 



m(/3;0)= lim 77; / vp{t;(j))dt, 

T->oo 1 Jq 

K{l3;(t>) ^ \mi sup |w/j(t; (/>)|, 

T^oo o<:t<T 



(43) 

(44) 
(45) 




a{[3;<t))= \im^\j^ I vp{t-(^Ydt, (46) 



(47) 



M4(/3;0) lim - / vp{t;(f))'^dt. 



In these conjectures, we always assume that the initial 
function (j){t) is taken such that the solution vp(t;(j)) is 
not a steady state solution, i.e., (j) satisfies the condition: 

- (j){t) + sin(27r/3(/i(t)) ^0, -1 < i < 0. (48) 

We then have the following conjectures. 

Conjecture 1 For any 4> G C([— 1,0],R) that satisfies 
(l48t. we have 



(49) 



Conjecture 4 For any (j) G C([— 1,0],M) satisfying pS)) . 
we have 



lim .,n ,N = 3. 
a^{l3; (j)) 

From Conjectures and [H the constant 
Xo= lim ^ 



(52) 



(53) 



is well defined for any </) e C([-l, 0], R) satisfying 
and independent of 4>. Conjecture [4] suggests that when 
/3 € / is sufficiently large, the density of the distribution 
of the time series vp{t; (j)) tends to a Gaussian with mean 
/i = 0, and standard deviation a{/3;(j)), but is truncated 
at ±K{(3;(I)). Therefore, let 

P^(z;0)- lim i / Hizail3;q^)-vp{t;cl)))dt (54) 

T->oo 1 Jq 

where H{-) is the Heaviside step function. Then Pp{z; (f) 
measures the probability that vp{t; </>) < za{l3; (p). 

Conjecture 5 Let he given by ([55)) and p{v; 0,1, Kq) 
be defined by (j39p . For any (f> S (7([— 1,0],R) satisfying 
(1151) . we have 



lim Pniz; 

3£l,l3->oo 



p(v; 0, 1, KQ)dv 



(55) 



for all z e 



for any {3 £ I . 



Conjecture[S]can be thought of as a Central Limit The- 
orem result for the irregular solutions of the differential 
delay equation 
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These conjectures were based on our numerical studies 
of the differential delay equation 



dv 
It 



= -v + F{v{t-l)) 



(56) 



with the nonlinear function F(v) taken to be a sinusoidal 
function. We suspect that the same results also hold for 
any bounded and oscillating nonlinear function such that 
the solution is "chaotic". For example, Figure [T^ shows 
the numerical results for the step function nonlinearity 



Fiv) 



H{sm{2TTPv)) - 



and the quasi-periodic nonlinearity 

F{v) = i sin(27r/3'u) + sin(2V27r^u) 



(57) 



(58) 



respectively. Therefore, the proposed conjectures may be 
universal for these deterministic "chaotic" dynamics. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper, we have studied a simple differential de- 
lay equation that displays a variety of behaviors, includ- 
ing chaotic solutions. 

In Section |ll] we carried out a complete bifurcation 
analysis for the steady state solutions. When ^ — 1, 
our analysis show that for any positive integer fc, there 
is an interval Ik = {ak,bk) that contains k + 1/4, such 
that when f3 £ Ik, there are two stable steady states 
(which are symmetric with respect to 0). Furthermore, 
{P,v) = (ak,^yk) are saddle node bifurcation points, 
and (/3,w) = (6fe,±Zfc) are Hopf bifurcation points. Ex- 
plicit expressions for ak,bk,yk, Zk are given in Section 
III Al When /3 increases past bk, two stable periodic so- 
lutions are generated at the Hopf bifurcation points. In 
additional to these regular solutions, when j3 > 0.85, the 
equation also has irregular solutions, which show chaotic 
behaviors. 

In Section lllli we numerically studied the probabilis- 
tic properties of the irregular (chaotic) solutions. Our 
simulations suggest that when f3 is large (7 = 1), the 
discrete sequences generated by irregular solutions 
v{t) (obtained by sampling each numerical solution every 
1000 steps, i.e., Vn = v{n x lOOOAi)) have the character 
of Gaussian distributed noise, but are truncated at the 
bound ±K which varies as ^-1/^ The variance of the 
time series also depends on (3 as (3~^^^. When (3 

is sufficiently large, the density of the distribution of the 
normalized solution approaches a quasi-Gaussian distri- 
bution ([M)) with parameter Kq = K/a ~ 5. 

In Section fllli the quasi-Gaussian distribution was ob- 
tained from the time series {w„} of solutions v{t) of Eq. 

with constant initial functions, and each solution is 
sampled every 1000 steps. We also noted that the sta- 
tionary density function is independent of the initial func- 
tion. We argue that the main results obtain in Section 
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FIG. 12. (Color online) (a) Upper bound K and standard 
deviation a for solutions of (|56|) with the step function (|57|) 
(hollow) and quasi-periodic function (solid) nonlinearities re- 
spectively. The dashed lines show the dependence K ~ /3~^/^ 
and o" ~ . In the case of the quasi-periodic function, 

the result for /3 = 4 (marked by an arrow) is exceptional be- 
cause the solutions are not chaotic (in fact, they are periodic 
solutions), (b) Normalized distributions obtained from the 
numerical solutions of H56p with a step function (blue circles) 
and quasi-periodic function (magenta triangles) nonlinearity 
respectively. The solid curve is the density function of the 
quasi-Gaussian distribution. Here /3 = 10. 



mil are independent of the sampling frequency. Thus, 
Figure [13] shows the bound K and standard deviation a, 
for different values of /?, of the time series when we 
sample the numerical solutions of Eq. every 1 step 
(f„ = v{n/S.t)). The results obviously show K ^ 
and a ^ fi~^^'^, as we have seen in Section [1111 Fig- 
ure 113b shows the distribution obtained from the time 
series obtained by sampling a solution with different fre- 
quencies (every 1 step, 500 steps, and 1000 steps, re- 
spectively). These simulation results suggested several 
conjectures (which we have been unable to prove) for the 
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probabilistic properties of the solutions of ([5]) as given in 
Section IVl 





FIG. 13. (Color online) (a) Upper bound K and standard de- 
viation a obtained from time series {wn.} obtained by sampling 
numerical solutions every 1 step, i.e., u„ = v{n/S.t). Dashed 
lines show the dependence K ~ and g ~ /3^^''^. Here 

7 = 1. (b) Normalized distributions obtained from time series 
obtained by sampling a solution with different frequencies (ev- 
ery 1 step (green circles), every 500 steps (black squares), and 
every 1000 steps (magenta triangles)). Here /3 = 10,7 = 1. 

Section HVl has shown that a Brownian motion like be- 
havior can be reproduced from the quasi-Gaussian dis- 
tributed solution of the differential delay equation. This 
deterministic Brownian motion shows behavior similar to 
that of experimentally observed Brownian motion, and 
therefore provides an alternative way to model appar- 
ently erratic behavior in nature. For example, the close 
to 50% efficiency exhibited in certain biological processes 
is very difficult to explain from a purely thermodynamic 
point of view [s^. The dynamical alternative presented 
here could afford another possible explanation, which 
could originate in a coherent, dynamical behavior at the 
molecular level of description. We feel that the applica- 



tion of the concept of a deterministic Brownian motion 
in modeling physical or biological phenomena that dis- 
play stochastic aspects will be of great interest in future 
studies. 

Finally, in Section |V] we have formulated five conjec- 
tures derived from our extensive numerical studies of this 
paper. We hope that these serve as a challenge to others. 

The significance of these results is, we feel, interesting. 
All experimental measurements typically exhibit fluctu- 
ations around some value, and it is customary (indeed 
the norm) to interpret these as "noise" and the implicit 
assumption is that these fluctuations are due to some ran- 
dom process that has no deterministic origin. The den- 
sity of the distribution of these fluctuations is, moreover, 
typically approximately Gaussian distributed but they 
are never truly Gaussian distributed (in the sense that 
the density is supported on the entire real line) but are 
always quasi-Gaussian in the sense that we have used it 
here. The numerical studies that we have presented lend 
strong circumstantial support to the alternative interpre- 
tation that what is typically held to be the signature of 
a random (i.e. non-deterministic) process could equally 
well be the signature of a completely deterministic pro- 
cess [Sll, [12]. The same implications were pointed out 
by Mackey and Tyran-Kamihska [l^ based on analytic 
computations in a similar situation. 



Appendix: Proof of Theorem [T] 

To prove theorem [Tl we first need the following two 
lemmas, which are obvious from the implicit function 
theorem, and the proofs are omitted. 

Lemma 2 Let (/JojWo) satisfy ([TT|) . and assume that 

1 - 27r/3ocos(27r/3ouo) 7^ 0. 

Then there is a function g{P), such that: 

(1) The function g{P) satisfies 5(/3o) — ^o? <^'^d 

g{(3) = sin(27r/3.g(/3)) 
for P in a neighborhood of /3q . 

(2) The function .g(/3) is differentiable, and 

dg[(i) 27rg(/3)cos(2^/3g(/3)) 



d/3 l-27r^cos(27r^g(^)' 

(3) The function g{l3) satisfies 
d{f3 cos(2^/35(/?))) cos(2^/3g(/3)) - 2^/3 



(A.l) 



dl3 1 - 27r/3cos(27r/3g(/3)) 

for (3 in a neighborhood of /3q. 
Lemma 3 Let (/3o,vo) satisfy (jlip . and assume that 

2Trvn cos(27r/3owo) 7^ 0. 
Then there is a function h{v), such that: 



(A.2) 
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(1) The function h{v) satisfies ft,(uo) = /3o> o.'rid 

h{v) — sin(27r/i(w)u) 
for V in a neighborhood of vq . 

(2) The function h{v) is differentiable, and 

dh{v) 1 — 2TTh{v) cos{2TTh{v)v) 
dv 2'KV cos{27rh{v)v) 

(3) The function h(v) satisfies 

d{h{v) cos{2Trh{v)v)) cos{2TTh{v)v) —2-Kh{v) 
dv 27ru cos(27r/i(u)w) 

for V in a neighborhood of vq . 



(A.3) 



(A.4) 



Proof of Theoreni[J\ Points (1), (2) and (3) are obvious. 

Point (4). Here, we only consider the positive solu- 
tions V, as the negative solutions are symmetric. The 
saddle node bifurcation points are given by solutions of 
the equations 



u sin(27r/3'y), 27r^ cos(27r^t;) = 1. 



(A.5) 



When w > is a solution of (jA.5|) . then sin(27r/3w) > 
and cos(27r/3u) > 0. Therefore solutions of (jA.Sp will 
always satisfy 

2fc7r < 27r/3u < 2fc7r + ^, fc G Nq. 



Hence, equations (jA.Sp arc equivalent to 

=1, 27r/3u = 2fc7r+arccos-^, fc G No. 
\27rp J 2ttp 

(A.6) 

Consequently, the bifurcation points are given by the so- 
lutions of 



2.fc = 2./3^1-(^^^ 

k 1 1 

V = — \ arccos 

/3 2ttP 27r/3 



arccos 



(A.7) 

From Lemma [1] these equations have a unique solution, 
which gives the saddle node bifurcation points {f3, v) = 
(afe,2/fc). 

Point (5). The Hopf bifurcation points arc given by 
the solutions of 



V — sin(27r/3u), 27r/3 cos(27r/3w) = secw. 



(A., 



Similar to our previous argument in Point (4), wc only 
consider positive solutions v that are given by the solu- 
tions of 




secw 



sec a; 

V ~ — I arccos 

13 27r/3 2ttI3 



27^/3 k e No, 
(A.9) 



which give the Hopf bifurcation points {bk, Zk). 
Point (6). Let 

F(/3,w) = v- sin(27r/3v). 

For any fc S N*, our previous arguments indicate that 
{ak,Vk) satisfies 



F(afc,yfe) = 0, 



and further 

dF{ak,yk) 
dp 



-27r?/fcC0s(27rafcZ/fc) 7^ 0. 



Thus, from Lemma[3l there is a function /3 = hk{v), such 
that ak = hk{yk), and it is differentiable in a neighbor- 
hood of Uk- 

Wc will show that the function hk{v) can be continued 
to the interval v g (0, yk]- If not, there is v* G (0, yk] and 
(3* such that 



F{p\v*) = 0, 



and 



dF{P*,v*) 
dp 



-27rw*cos(27r/3*w*) = 0, 



which implies cos(27r/3*i;*) = 0. Therefore we should 
have 

{v*f = (sin(27r/3*w*))2 = 1. 

However, this is impossible since yk < 1- Thus, wc con- 
clude that the function hk{v) can be continued to the en- 
tire interval {0,yk], and further that cos{2TThk{v)v) > 
for any v S (0,yfc]. 

Next, we will show that hk{v) is a decreasing function 
for V € (0, yk]- From (jA.4|) . we have 



< 0. 



d{hk{v) cos{2Trhk{v)v)) cos{2nhk{v)v) — 2'Khk{v) 
dv 2t:v cos{2T:hk{v)v) 

Thus, for V € {0,yk), wc have 

2Trhk{v) cos{2T:hk{v)v) > 2TrakCos(2'Kakyk) = 1- 



Note that hk{v) > oo = 1/{2tt) and cos{2TThk{v)v) > 0. 
Therefore, from Lemma [3] 



dhk{v) 1 — 2iThk{v) cos{2Trhk{v)v) 



dv 



2'nv cos{2'Khk {v)v) 



< 0. 



Now, the function (3 = hk{v) is well defined and de- 
creasing for V € (0, i/fc]. Thus, the inverse function, de- 
noted by w = fkiP), is also well defined, continuous at 
(3 e [afc,oo), and such that v = ±fk{(3) satisfy JTT]) and 
(P^. From (22), it is easy to conclude that the steady 
state solutions v{t) = ztfk{P) are unstable, and Point (6) 
is proved. 
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Point (7). For any k € No, the Hopf bifurcation point 
{bk, Zk) satisfies 

1 — 2TTbk cos(27r6fcZfc) = 1 — sec w > 0. 

Therefore, we can apply Lemnia[21 and there is a function 
gkiP) such that v = gk{P) satisfies (HI]), and Zk — gk{bk)- 
When f3 > bk (> ao), from (jA.2|) in Lemma [21 we have 

d(l-2^/3cos(2^/3gfc(/3)))^ 
df3 

= -47r(cos(27r/3gfc(/3)) - 27r/3) > 0. 

Thus, we have dF{f3, v)/dv ^ for (3 > bk and the func- 
tion gfc(/3) can be continued to /3 S {bk, oo), and the 
steady state solutions v{t) = ±gk{l3) are unstable. 

When f3 < bk, we will show that the function gk{P) can 
be continued to /3 € (a^, bk). If not, there is /3* G {ak,bk) 
such that 

F{/3*,v*) = 

and 

^^'^^l*^'"*^ = 1 _ 27r/3* cos(27r^*w*) = 0. 
ov 

Therefore, we must have /3* — ak' for some k' e Nq. 
However, ak is the maximum of such values that are less 
than bk, and thus we must have /3* = ak. Therefore, 
the function gkiP) can be continued to /3 G {ak,bk), and 
gk{cLk) = Vk- These arguments show that the function 
gk{P) is well defined in the interval [ak, bk), and satisfies 
Vk = gk{ak), Zk = gk{bk). 

Now, we only need to show that when ak < 13 < bk, 
the steady state solutions v{t) = ±gk{P) are stable. Since 



(1 — 27r/3 cos(27r/3(7fe(/3)))^ is increasing with respect to (3, 
and 

f 1 - 27rafe cos(2nakgk{a,k)) = 0, 

1^ 1 — 27r5fe cos(27r6fegfe(6fe)) = 1 — secw > 0, 

we have 

< 1 - 2nl3 cos{2tt 13 gkil3)) < 1 - secw, 

i.e., 

secw < 27r/3cos(27r/35fe(/3)) < 1 

for any /3 E {ak, bk). Therefore the steady state solutions 
v{t) = gk{P) are stable, and Point (7) is proved. 

Point (8) is obvious from the above arguments, and the 
theorem is proved. 
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